#ifndef __ROUTINE_NUMERIC__
#define __ROUTINE_NUMERIC__

#include "stdio.h"
#include "stdlib.h"
#include "complex.h"
#include "math.h"

#include "fftw3.h"

#define c_vacc (137.036)       // in a.u. units
#define PI (3.1415926535897932)

double det(double *A, int n);

/**
 * Doing integrations
 * 
 */
double complex Intl_grid_fun_1(double complex *grid_fun,
			       int Nx, int Ny, int Nz,
			       double hx, double hy, double hz);

void Intl_grid_coef(double *coef, int Nx, int Ny, int Nz);


/**
 *   1. FFTW Code for -\laplacian u= 4*pi*rho
 *   
 */
void poissonFFT(double complex *v_H, double complex *rho, 
		int Nx, int Ny, int Nz, double hx, double hy, double hz);

/**
 * 2. Second order gradients for grid functions
 * Using periodic condition 
 */

//template<class typename>
void gradient_real_2(double *f, double *fx, double *fy, double *fz, 
		     int Nx, int Ny, int Nz, double hx, double hy, double hz);

void gradient_cplx_2(double complex *f, 
		     double complex *fx, double complex *fy, double complex *fz, 
		     int Nx, int Ny, int Nz, double hx, double hy, double hz);

/**
 * 
 * 
 */
void gradient_of_complex_function(double complex *f, 
				  double complex *fx, 
				  double complex *fy, 
				  double complex *fz, 
				  int Nx, int Ny, int Nz, 
				  double hx, double hy, double hz);

void gradient_of_real_function(double  *f, double  *fx, double  *fy, double  *fz, 
			       int Nx, int Ny, int Nz, 
			       double hx, double hy, double hz);

void gradient_of_wfs(double *wfs, double *gradient_wfs, int st_num, 
		     int Nx, int Ny, int Nz, 
		     double hx, double hy, double hz);

void gradient_hessian_eta_xc(double complex *eta_xc,
			     double complex *gradient_eta_xc, 
			     double complex *hessian_eta_xc, 
			     int Nx, int Ny, int Nz, 
			     double hx, double hy, double hz);

/* /\** */
/*  * 3. Fourth order */
/*  * Using periodic condition  */
/*  *\/ */

/* template<class typename> */
/* void grad_of_gridfun_4(typename *f, typename *fx, typename *fy, typename *fz,  */
/* 		       int Nx, int Ny, int Nz, double hx, double hy, double hz); */

/**
 * interpolate the field A_EM from FE mesh to molecule grid
 * 
 */

void mesh_2_grid(double complex *AEM_1, double complex *AEM_2,  
		 double complex *AEM_3,   
		 double *x1, double *y1, double *z1, int Nx, int Ny, int Nz,
		 double *sampleX, double *sampleY, double *sampleZ, int N_sample, 
		 double complex *A_EM, 
		 double *Node, int *Element, int *Element_Edge, 
		 double *Length_of_Edge, 
		 int *PointInTet);

/**
 * interpolate the current density j to FE mesh
 * 
 */

void grid_2_mesh(double complex *del_j ,    // out put FEM function
		 double *Node, int *Element, int *Element_Edge, 
		 int N_edge, double *Length_of_Edge,
		 double complex *del_j_1,   // input variables
		 double complex *del_j_2, 
		 double complex *del_j_3,
		 double *x1, double *y1, double *z1, int Nx, int Ny, int Nz,
		 double *sampleX, double *sampleY, double *sampleZ, int N_sample, 
		 int *NodeInSamples);


#endif
